#Run transformation2.R first

#Working directory to file directory

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))


#Create Plots folder, if it does not exist already
dir.create(file.path(getwd(), "Plots"))

#Install Required Packages


packages <- c("broom", "sjPlot", "estimatr", "tidyverse", "ggplot2", "patchwork", "stargazer", "svglite", "ggeffects")

# Check if each package is installed, install if missing, and load it
for (pkg in packages) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg)
    require(pkg, character.only = TRUE)
  }
}



#####
#Table 2
#####
Tab2M <- lm(checkB ~ officeF, data = data2)
Tab <- ggpredict(try, terms = "officeF")

#####
#Figure 3
#####


#violation
violation  <- lm_robust(choice ~ violation, 
                       data = data2, clusters = caseid)

#tab_model(violation, pred.labels = c("Compliance", "Norm Violation"),
#          dv.labels = "Defection Rates",
#          string.ci = "Conf. Int (95%)",
#          string.p = "P-Value")

predictions <- predict(violation, newdata = data2, se.fit = TRUE)
data2$violation_fit <- predictions$fit
data2$violation_se.fit <- predictions$se.fit
data2$violation_lower_ci <- data2$violation_fit - 1.96 * data2$violation_se.fit
data2$violation_upper_ci <- data2$violation_fit + 1.96 * data2$violation_se.fit

plot_data2 <- data2 %>%
  group_by(violation) %>%
  summarise(
    Mean_Prediction = mean(violation_fit),
    Lower_CI = mean(violation_lower_ci),
    Upper_CI = mean(violation_upper_ci),
    .groups = 'drop'
  )

violationP <- ggplot(plot_data2, aes(violation, Mean_Prediction, fill = factor(violation))) + geom_bar(stat="identity", position=position_dodge(), alpha = .7) + geom_point(size = 2) +
  geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI), width = .1, size = 1) + xlab("") + ylab("Defection Rates") +
  theme_minimal() +
  scale_x_continuous(breaks = c(0, 1), labels = c("Compliance", "Norm Violation")) +
  theme(legend.position = "none",
        text = element_text(family = "serif", size = 12)) + scale_y_continuous(breaks = c(0, .1, .2, .3, .4), 
                                                                               labels = c("0 %", "10 %", "20 %", "30 %", "40 %"), limits = c(0, .5)) +
  coord_flip() + scale_fill_manual(values = c("0" = "lightblue3", "1" = "cornflowerblue"))


norm_type <- lm_robust(choice ~ normF + violation + normF*violation, 
                       data = data2, clusters = caseid)


pred_data <- expand.grid(
  normF = unique(na.omit(data2$normF)),
  violation = unique(na.omit(data2$violation))
)

pred_data$Predictions <- predict(norm_type, newdata = pred_data, se.fit = TRUE)$fit
pred_data$SE <- predict(norm_type, newdata = pred_data, se.fit = TRUE)$se.fit
pred_data$Upper <- pred_data$Predictions + 1.96 * pred_data$SE
pred_data$Lower <- pred_data$Predictions - 1.96 * pred_data$SE


normtype <- ggplot(pred_data, aes(x = normF, y = Predictions, fill = factor(violation))) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7, alpha = 0.7) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = .1, size = 1, position = position_dodge(width = 0.8)) +
  geom_point(position = position_dodge(width = 0.8), size = 2) +
  labs(x = "", y = "Defection Rates", title = "") +
  theme_minimal() + theme(text = element_text(family = "serif", size = 12),
                          legend.position = "null") + 
  scale_fill_manual(values = c("0" = "lightblue3", "1" = "cornflowerblue")) +
  scale_y_continuous(breaks = c(0, .1, .2, .3, .4), 
                     labels = c("0 %", "10 %", "20 %", "30 %", "40 %"), limits = c(0, .5)) +
  geom_text(data = pred_data, 
            aes(label = ifelse(violation == 0, "", "Violation"), y = 0), 
            position = position_dodge(width = 0.8), 
            family = "serif", color = "black", size = 3.2, hjust = -0.2,
            fontface = "bold") + coord_flip()



office <- lm_robust(choice ~ violation + officeF + officeF*violation, data = data2, clusters = caseid)

pred_data <- expand.grid(
  officeF = unique(na.omit(data2$officeF)),
  violation = unique(na.omit(data2$violation))
)

pred_data$Predictions <- predict(office, newdata = pred_data, se.fit = TRUE)$fit
pred_data$SE <- predict(office, newdata = pred_data, se.fit = TRUE)$se.fit
pred_data$Upper <- pred_data$Predictions + 1.96 * pred_data$SE
pred_data$Lower <- pred_data$Predictions - 1.96 * pred_data$SE

pred_data$officeF <- factor(pred_data$officeF, levels = c("Local", "State", "Federal"))


officeP <- ggplot(pred_data, aes(x = violation, y = Predictions, fill = officeF)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7, alpha = 0.7) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = .1, size = 1, position = position_dodge(width = 0.8)) +
  geom_point(position = position_dodge(width = 0.8), size = 2) +
  labs(x = "", y = "Defection Rates", title = "") +
  theme_minimal() + theme(text = element_text(family = "serif", size = 12),
                          legend.position = "null") + 
  scale_fill_manual(values = c("Local" = "lightblue3", "State" = "cornflowerblue", "Federal" = "skyblue2")) +
  scale_y_continuous(breaks = c(0, .1, .2, .3, .4), 
                     labels = c("0 %", "10 %", "20 %", "30 %", "40 %"), limits = c(0, .5)) +
  scale_x_continuous(breaks = c(0, 1), labels = c("Compliance", "Norm Violation")) +
  geom_text(data = subset(pred_data, violation != 0), 
            aes(label = officeF, y = 0), 
            position = position_dodge(width = 0.8), 
            family = "serif", color = "black", size = 3.2, hjust = -0.2,
            fontface = "bold") + coord_flip()




#factual manipulation check

data2 = data2 |> mutate(checkB = case_when(Check_flag == "Pass" ~ 1,
                                           Check_flag == "Fail" ~ 0,
                                           TRUE ~ NA_real_))




#Figure 3

violationP <- violationP + theme(text = element_text(size = 15))
normtype <- normtype + theme(text = element_text(size = 15))
officeP <- officeP + theme(text = element_text(size = 15))


# Combine the plots with the specified layout
combined_plot <- (violationP / officeP) +
  plot_layout(heights = c(1, 2)) | normtype


svg("plots/Figure3.svg", height = 6, width = 9)
combined_plot + plot_annotation(title = expression(italic("Including Policy Tradeoff"))) &
  theme(plot.title = element_text(hjust = 0.5))
dev.off()


